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Introduction 

The unsteady evolution of three-dimensional synthetic jet into quiescent air is studied by 
time-accurate numerical simulations using a second-order accurate mixed explicit- 
implicit fractional step scheme on Cartesian grids. Both two-dimensional and three- 
dimensional calculations of synthetic jet are carried out at a Reynolds number (based on 
average velocity during the discharge phase of the cycle Vj, and jet width d) of 750 and 
Stokes number of 17.02. The results obtained are assessed against PIV and hotwire 
measurements provided for the NASA LaRC workshop on CFD validation of synthetic 
jets. 

Numerical Methodology 

The evolution of zero-net mass-flux synthetic jet from a cavity into quiescent air is 
modeled by the unsteady, incompressible Navier-Stokes equations, written in tensor form 
as 
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where the indices, i =1, 2, 3, represent the x, y and z directions, respectively; while the 
velocity components are denoted by u (u 0, v (M 2 ) and w (M 3 ), respectively. The equations 
are nondimensionalized with the appropriate length and velocity scales where Re 
represents the Reynolds number. The Navier-Stokes equations are discretized using a 
cell-centered, collocated (non-staggered) arrangement of the primitive variables (u ,p). In 

addition to the cell-center velocities (m ), the face-center velocities,!/ , are also computed. 
Similar to a fully staggered arrangement, only the component normal to the cell-face is 
calculated and stored. The face-center velocity is used for computing the volume flux 
from each cell. The advantage of separately computing the face-center velocities has been 
discussed in the context of the current method in Ye et al. [1]. The equations are 
integrated in time using the fractional step method. In the first step, the momentum 
equations without the pressure gradient terms are first advanced in time. In the second 
step, the pressure field is computed by solving a Poisson equation. A second-order 
Adams-Bashforth scheme is employed for the convective terms while the diffusion terms 
are discretized using an implicit Crank-Nicolson scheme which eliminates the viscous 
stability constraint. The pressure Poisson equation is solved with a Krylov-based 
approach. 
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The solver uses a multi-dimensional ghost-cell methodology to incorporate the effect 
of the immersed boundary on the flow; however, the absence of curvilinear topology in 
the geometry makes this approach unnecessary in the current study. Care has been taken 
to ensure that the discretized equations satisfy local and global mass conservation 
constraints as well as pressure-velocity compatibility relations. The code has been 
rigorously validated by comparisons against established experimental and computational 
data. Details have been presented elsewhere [2, 3]. 

Implementation and Case Specific Details 

The present study models the flow inside the cavity using a pulsatile velocity boundary 
condition, v = v a sin(ftl) prescribed at the bottom of the cavity (see Fig. 1) in order to 
generate a natural flow at the slot exit. The shape of the cavity is approximated to be a 
rectangular box without taking into consideration the finer details that make up the 
interior of the cavity. However, the geometrical and flow parameters used in the current 
study are chosen based on a scaling analysis of various parameters involved like slot 
width (d), slot height (h), cavity width (W), cavity height (H), the diaphragm vibrating 
frequency (f), etc.. For instance, the slot size is chosen such that the ratio h/d and W/d 
match those used in the experiments; however, H/d of 4.95 used in the calculations is not 
matched in the experiments. The Reynolds number (Re) in this work is defined based on 
average jet velocity during the discharge phase of the cycle (Vj) and jet width 
(d)i.e.. Re = Vjd / v . Average jet velocity (Vj) in the numerical calculations is set equal to 

the value of 10.5 m/s, obtained by averaging the velocity provided by the LDV 
measurements during the discharge phase of the cycle. Because of the uncertainty in Vj 
reported by different measurements, a nominal Re of 750 corresponding to the lower 
bound of Vj calculated from the measurements is chosen in the present investigation. 
Experimental values of ft) = 2794 rad/s, d = 1.27 mm and V = 1.5527xl0' 5 m/s are 

matched to give a Stokes number ( S = V cod 2 / V ) of 17.02 in the computations. Various 
cases considered in the present study and the corresponding flow parameters are detailed 
in Table 1. 
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3D 

17.02 
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0.00116144 

2. 

750 

2D 

17.02 

30d x 30d 

132x220 

14,000 

0.00116144 


Table l.Various cases considered in the study and their flow parameters. 


Figure 1 shows the schematic of the computational domain and the boundary conditions 
used in the computations. An outflow velocity boundary condition is prescribed at the 
left, right and top boundaries that allow them to respond freely to the flow created by the 
jet. In 3D calculations, periodic boundary conditions are prescribed in the span wise (z) 
direction. Figure 2 shows an x-y slice of a typical 3D mesh used in the region near the 
slot in the computations. Grids used in the current work are non-uniform in both x- and y- 
directions, and uniform in the span wise (z) direction in the case of 3D calculations. 
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Sufficient clustering is provided in the slot-region along x- and y-dircctions to resolve the 
vortex structures that form at the slot exit, as well as the shear layer in the slot. Typically 
32 grid points clustered using a cosine-hyperbolic distribution are used across the slot. In 
3D calculations, the three-dimensionality in the solution is instigated by introducing a 
small sinusoidal spatial perturbation in the z-component of velocity (w) over a few 
hundred time- steps in the first cycle. For all cases presented here, the first two cycles of 
calculations are not included in the averaging process to eliminate transient effects, and 
the next two to three cycles are used in the process. 

Results 

Figure 3 shows the isosurfaces of vorticity magnitude obtained from three-dimensional 
calculation before the onset of full three-dimensionality in the solution. It is clear from 
the figure that the flow is dominated by counter-rotating vortex pairs. Figure 4 shows the 
isosurfaces of vorticity magnitude for the same solution after the onset of full three- 
dimensionality. Figure 4 shows two pairs of rib vortices along the span obtained from the 
3D simulation. The plot of phase-averaged v-velocity component vs. phase angle <j) at the 
point (. x. , y) = (0 mm, 0.1 mm) shown in Figure 5 describes the procedure involved in 
aligning the CFD data with PIV data. Phase angle at which V avg = (V max + V m m)/2 intersects 
the curve is made 340° by applying the required phase shift. Figure 6 shows the plot of 
phase-averaged v-velocity component vs. phase angle (f) at the point (x, y) = (0 mm, 0.1 
mm) for the two calculations and it can be seen that the CFD data is reasonably aligned 
with the PIV data. Plot of time-averaged u- and v-velocity components at the point (x, y) 
= ( 0 mm, 0.1 mm) across the slot region is shown in Figure 7. Figure 8 and 9 show the 
comparison of phase-averaged u- and v-velocity components at the same point at (/) = 90° 
and <j> = 270° respectively. Three-dimensional calculations at Re > 1000 corresponding 
to the upper bound of V, reported in the experiments are being carried out and since the 
solutions at these Reynolds numbers were not converged at the time of preparation of this 
report, they would be presented in detail at the workshop. 
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Figure 1: Schematic of the computational 
domain and boundary conditions. 



Figure 2: A typical grid in the slot region 
used in the computations. 



Figure 3: Isosurfaces of vorticity magnitude Figure 4: Isosurfaces of vorticity magnitude 
2 nd cycle, Re = 750, S = 17.02. 3 rd cycle, Re = 750, S = 17.02 



Figure 5: Plot of phase-averaged v vs. phase Figure 6: Plot of phase-averaged v vs. phase 

angle before aligning CFD data with PIV data. angle after aligning CFD data with PIV data. 
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Figure 7 : Plot of time-averaged u- and v-velocities along the horizontal line 
y = 0.7874 (~lmm in the experiments). 



Figure 8: Plot of phase-averaged u- and v-velocities along the horizontal line 
y = 0.7874 (~lmm in the experiments) at (f) = 90° degrees. 



Figure 9: Plot of phase-averaged u- and v-velocities along the horizontal line 
y = 0.7874 (~lmm in the experiments) at (f>- 270° degrees. 
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